######## Attitudes to crime and punishment in England and Wales, 1964-2023
######## 003. Analysis_Responsiveness

#### Set up

library(tidyverse)
library(gt)

read.csv("/Users/matteo/Desktop/Prison Polling/LatentTrendsAllAdults.csv") %>%
  full_join(., read_csv("/Users/matteo/Desktop/Prison Polling/CrimeRate.csv", col_select = c(1:8))) %>%
  arrange(Year) %>%
  mutate(across(starts_with("CSEW"), ~ . /100),
         across(starts_with("Police"), ~ (. /Population)),
         HomicidePerMill = HomicidePerMill/1000000,
         Adjusted_HomicidePerMill = Adjusted_HomicidePerMill/1000000) %>%
  select(-c(Population,X)) -> Trends


# Figure 6: Victimisation rates in England and Wales, 1960-2023
dropLeadingZero <- function(l){stringr::str_replace(l, '0(?=.)', '')}
options(scipen = 6)
Trends %>%
  select(c(1,6:10)) %>%
  pivot_longer(-Year) %>%
  mutate(name = case_match(name, 'HomicidePerMill'~'Murder rate (police records)','Police_Property'~'Property crime (police records)',
                           'Police_VATP'~'Violent crime (police records)','CSEW_Property'~'Property crime (CSEW)',
                           'CSEW_VATP'~'Violent crime (CSEW)')) %>%
  mutate(name = factor(name, levels = c('Violent crime (CSEW)','Property crime (CSEW)','Violent crime (police records)','Property crime (police records)',
                                        'Murder rate (police records)'))) %>%
  drop_na(value) %>%
  ggplot(aes(x = Year, y = value, group = name)) +
  geom_line() + facet_wrap(facets = vars(name), nrow = 3, ncol =2 , scales = "free_y") +
  xlab(NULL) + ylab(NULL) +
  MyThemes::theme_base() +
  scale_y_continuous(labels = dropLeadingZero) +
  theme(plot.title = element_text(face = 'bold', hjust = 0.5))



###### Table 4: Bivariate correlations between crime and attitudes

corr_analysis <- function (var1, var2){
  cf <- cor.test(Trends[[var1]], Trends[[var2]])$estimate[[1]]
  pv <- cor.test(Trends[[var1]], Trends[[var2]])$p.value[[1]]
  n <- sum(complete.cases(bind_cols(Trends[[var1]],Trends[[var2]])))
  return(c(cf,pv,n,var1,var2))
}

lapply(c("CrimeConcern","Punitiveness","MIP","DeathPen1"), function(var1) {
  lapply(c("HomicidePerMill","Police_Property",
           "Police_VATP","CSEW_Property","CSEW_VATP"),
         function(var2){corr_analysis(var1, var2)})}) -> results

data.frame(matrix(unlist(results), nrow=20, byrow=TRUE),stringsAsFactors=FALSE) %>%
  rename(r = X1, p = X2, n = X3, var1 = X4, var2 = X5) %>%
  mutate(r = as.numeric(r), p = as.numeric(p)) %>%
  mutate(var1 = case_match(var1, 'CrimeConcern' ~ 'Fear', 'Punitiveness' ~ 'Punitiveness', 'MIP' ~ 'Prioritisation','DeathPen1' ~"Death Penalty"),
         var2 = case_match(var2, 'HomicidePerMill'~'Murder rate (recorded)','Police_Property'~'Property crime (recorded)',
                           'Police_VATP'~'Violent crime (recorded)','CSEW_Property'~'Property crime (survey)',
                           'CSEW_VATP'~'Violent crime (survey)')) %>%
  pivot_wider(values_from=c(r,p,n),names_from=var1) %>%
  select(var2,r_Fear,n_Fear,
         r_Prioritisation,n_Prioritisation,
         r_Punitiveness,n_Punitiveness,
         #`r_Death Penalty`,`p_Death Penalty`,`n_Death Penalty`
  ) %>%
  gt(rowname_col = "var2") %>%
  tab_spanner_delim(delim = "_", reverse = TRUE) %>%
  tab_header(title = md("**Table 4: Bivariate correlations between crime and attitudes**")) %>%
  fmt_number(decimals = 2) %>%
  tab_style(style = cell_text(style = "italic"), locations = cells_column_labels()) %>%
  cols_align(align = "center", columns = c(2:7)) %>%
  fmt_number(pattern = "{x} ***", columns = c(2), rows = c(1,2,4,5)) %>%
  fmt_number(pattern = "{x} ***", columns = c(4), rows = c(1,2,5)) %>%
  fmt_number(pattern = "{x} **", columns = c(6), rows = c(1)) %>%
  fmt_number(pattern = "{x} ***", columns = c(6), rows = c(2:5)) %>%
  cols_width(1 ~ pct(30), 2:7~ pct(70/6)) %>%
  tab_source_note(source_note = "Notes: Pearson product moment correlations.")


# Appendix 13: Alternative homicide counts

lapply(c("CrimeConcern","Punitiveness","MIP","DeathPen1"), function(var1) {
  lapply(c("Adjusted_HomicidePerMill"),
         function(var2){corr_analysis(var1, var2)})}) |> print()

rm(corr_analysis, results)



###### Table 5a: Maximum-likelihood factor analysis: 

Trends %>% mutate(across(!(Year), ~ scale(.))) -> lmready

efa_1 <- factanal(x = drop_na(lmready[,c(6,7,9,10)]),
                  factors = 1,
                  scores = "regression")
efa_1$loadings

tibble(Index = efa_1$scores, 
       Year = lmready$Year[!(is.na(lmready$CSEW_Property) | is.na(lmready$Police_Property))]) |>
  ggplot(aes(x = Year, y = Index)) + geom_line()


##### Table 5b: Bivariate regressions between crime and attitudes

lmready <- full_join(lmready, 
                 tibble(Index = efa_1$scores, 
                        Year = lmready$Year[!(is.na(lmready$CSEW_Property) | is.na(lmready$Police_Property))]),
                 by = 'Year')
broom::tidy(lm(CrimeConcern ~ Index, data = lmready))
broom::tidy(lm(Punitiveness ~ Index, data = lmready))
broom::tidy(lm(MIP ~ Index, data = lmready))



# Appendix 13: Alternative homicide counts

read.csv("/Users/matteo/Desktop/Prison Polling/LatentTrendsAllAdults.csv") %>%
  full_join(., read_csv("/Users/matteo/Desktop/Prison Polling/CrimeRate.csv", col_select = c(1:8))) %>%
  arrange(Year) %>%
  mutate(across(starts_with("CSEW"), ~ . /100),
         across(starts_with("Police"), ~ (. /Population)),
         HomicidePerMill = HomicidePerMill/1000000,
         Adjusted_HomicidePerMill = Adjusted_HomicidePerMill/1000000) %>%
  select(-c(Population,X)) %>% 
  mutate(across(!(Year), ~ scale(.))) -> lmready

efa_1 <- factanal(x = drop_na(lmready[,c(11,7,9,10)]),
                  factors = 1,
                  scores = "regression")
efa_1$loadings

tibble(Index = efa_1$scores, 
       Year = lmready$Year[!(is.na(lmready$CSEW_Property) | is.na(lmready$Police_Property))]) |>
  ggplot(aes(x = Year, y = Index)) + geom_line()

lmready <- full_join(lmready, 
                     tibble(Index = efa_1$scores, 
                            Year = lmready$Year[!(is.na(lmready$CSEW_Property) | is.na(lmready$Police_Property))]),
                     by = 'Year')
broom::tidy(lm(CrimeConcern ~ Index, data = lmready))
broom::tidy(lm(Punitiveness ~ Index, data = lmready))
broom::tidy(lm(MIP ~ Index, data = lmready))
